Divergence of the Thermal Conductivity in Uniaxially Strained Graphene 
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We investigate the effect of strain and isotopic disorder on thermal transport in suspended 
graphene by equilibrium molecular dynamics simulations. We show that the thermal conductiv- 
ity of unstrained graphene, calculated from the fluctuations of the heat current at equilibrium is 
finite and converges with size at finite temperature. In contrast, the thermal conductivity of strained 
graphene diverges logarithmically with the size of the models, when strain exceeds a relatively large 
threshold value of 2%. An analysis of phonon populations and lifetimes explains the divergence of 
the thermal conductivity as a consequence of changes in the occupation of low-frequency out-of-plane 
phonons and an increase in their lifetimes due to strain. 

PACS numbers: 65.80.Ck, 63.22.Rc, 05.60.Cd 
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I. INTRODUCTION 

The combination of light weight, strong covalent 
bonds, and low dimensionality gives carbon nanostruc- 
tures, such as graphene and nanotubes, superior me- 
chanical and thermal properties^, making them interest- 
ing candidate materials for thermal management^ and 
phononics applications^Ti^. Extremely high, possibly di- 
vergent, thermal conductivity (k) of a two-dimensional 
(2D) phonon gas was predicted by Klemens and Pe- 
draza ten years before single layer graphene was iso- 
lated for the first time^. The first measurements in sus- 
pended graphene at about room temperature confirmed 
Klemens and Pedraza's predictions, reporting values of 
K in the range 3000 - 5800 Wm^^K-^i^i. Later ex- 
periments found K w 2500 Wm~^K~^at 350 K, and 
« 1400 Wm"^K"iat 500 K^. Recent measurements 
on suspended graphene in vacuum yield k in the range 
2600 - 3100 Wm-^K-^at 350 K^. 

In spite of the efforts to refine these measurements, 
an accurate determination of k remains a tough exper- 
imental challenge^. A fundamental reason for such dif- 
ficulty is that heat transport in graphene is very sen- 
sitive to defects and experimental conditions. For ex- 
ample, when graphene is supported on a substrate k is 
reduced to w 600 Wm-^K-^at room temperature due to 
phonon scattering from the substrate^^. k can vary by 
as much as 50% at room temperature as a function of 
the isotopic composition of graphenoii, and was found 
to be sensitive also to the lateral dimension of measured 
patches^ and to the presence of wrinkles, which may 
lower K by 30%^^. This high sensitivity can be advan- 
tageous because it offers the possibility to manipulate n 
in graphene-based devices either by tuning the concen- 
tration of defects and mass disorder during growth or 
by imposing controlled external conditions. High and 
tunable thermal conductivity for a single sheet of atoms 
opens up the possibility of application in a range of ther- 
mal management devices, from high-power electronics^^ 
all the way to thermoelectric applications^^. 

In pristine graphene at room temperature, heat is con- 
ducted almost exclusively by phonons^, so we can fo- 



cus on lattice thermal conductivity. Both lattice dynam- 
ics (LD) calculationai^"— and molecular dynamics (MD) 
studie o^^'^"" — indicate that k of suspended graphene con- 
verges with system size, in contrast with ideal 2D mod- 
els, for which k diverges logarithmicall}^"— . Flexural 
phonons (out-of-plane vibrational modes) play a decisive 
role both as heat carriers and as scatterers. Neverthe- 
less, some of these works report very diverse numerical 
results, stemming from the use of different methods rely- 
ing on different approximations, and from the choice of 
various interatomic potentials. Theoretical studies sug- 
gest that it is possible to control k in graphene by ap- 
plying mechanical (tensile) strain, however simulations 
have given contradictory results. Recent ab initio LD 
calculations showed that, while k is finite for unstrained 
suspended graphene, it diverges when tensile strain is 
appliedi^. In contrast, MD results point in the oppo- 
site direction, indicating a reduction of k upon strain2£. 
Even though discrepancies between LD and MD results 
are expected, as the two methods rely upon different ap- 
proximations2L^, it is unusual to get such differences 
in trends. In fact, in LD calculations anharmonic inter- 
actions are usually truncated at the first order, while in 
MD simulations quantum effects cannot be taken into ac- 
count, so phonon populations obey to classical statistics. 
Both approximations conspire to make kld larger than 

HMD- 

In this work we report the results of equilibrium molec- 
ular dynamics (EMD) simulations of heat transport in 
suspended graphene as a function of strain and isotopic 
mass disorder. We begin by investigating size conver- 
gence of K in isotopically pure unstrained graphene, and 
then we study the effects of mechanical strain and iso- 
topic disorder. Our goal is to verify whether the di- 
vergence, predicted by LD calculations on isotropically 
strained graphene, also occurs upon uniaxial strain at fi- 
nite temperature. We also probe how the combination of 
strain and mass disorder affects k. A microscopic inter- 
pretation of the results is provided in terms of phonon 
populations and lifetimes computed at finite tempera- 
ture. 
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II. METHODS 

We compute the thermal conductivity of graphene by 
EMD simulations in models with periodic boundary con- 
ditions. We use the TersofF empirical potential^^ recently 
re-parametrized to accurately reproduce the vibrational 
properties of carbon nanostructuresi^. Anharmonic LD 
calculations employing this set of parameters result in a 
thermal conductivity k = 3500 Wm^^K^^for a 10 /im 
graphene flake at T = 300 K, well within the range of 
experimental measurements. MD production runs are 
performed in the microcanonical ensemble, starting from 
initial configurations equilibrated at the target tempera- 
ture^. Temperatures between 300 and 1000 K are con- 
sidered. The equations of motion are integrated with a 
1 fs time step. The cell parameters are optimized at the 
simulation temperature to achieve zero stress in the (zig- 
zag) direction, perpendicular to the strained one. 

Following linear response theory, n is computed from 
the integral of the autocorrelation function of the heat 
flux J(t) in a microcanonical simulation, according to 
the Green-Kubo formulap^i*^ 

'^"/^ -T^I J™ Ti I {Mt')MO))dt', (1) 

Kb-L t-fooV^oo V Jq 

where fcs is Boltzmann's constant, T is the temperature 
and V the volume, which is here defined as the surface 
area of the graphene foil times a nominal thickness of 
3.35 A. In practice, k is taken as the stationary value 
of Eq. ^ before it drifts due to accumulated statistical 
noise. Although k is in general a tensor, the hexagonal 
symmetry of graphene yields k^x = i<-yy = and K^y = 0. 
The limits to infinite time and infinite volume in Eq. ([T]) 
indicate that size convergence and phase space sampling 
have to be carefully considered. This aspect is particu- 
larly important for low-dimensional systems, for which 
transport coefficients usually diverge^^"— i^. Therefore, 
investigating size and time convergence is not merely a 
technical aspect, but it brings important physical insight. 
In order to effectively sample the phase space and achieve 
statistical accuracy in evaluating Eq. ([1]), each reported 
value of K is obtained by averaging over at least 20 inde- 
pendent simulations of at least 60 ns. 



III. RESULTS AND DISCUSSION 

A. Thermal conductivity of unstrained suspended 
graphene 

To check size convergence we perform simulations with 
approximately square supercells of increasing size. We 
consider systems made of from 240 to 3 • 10^ atoms. 
The smallest supercell is 26 x 25 and the largest one 
880 X 877 A^ . Fig. [T]shows the calculated thermal conduc- 
tivity as a function of the number of atoms in the simu- 
lation cell. The anisotropy between the in-plane thermal 
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FIG. 1. (Color online) Thermal conductivity of graphene as a 
function of the number of atoms in the simulation cell. Con- 
verged value K = 1015 ± 120 Wm~^K~^is achieved for a cell 
with 17200 atoms. Inset: time decay of the normalized heat 
flux autocorrelation functions (HFACF) for 4400, 17200, and 
287232 atom cells. AU HFACF decay faster than 1/t. 



conductivities seen for the smallest cell (240 atoms) is 
due to an uneven and insufficient sampling of the vibra- 
tional modes in the two directions. As the cell size is in- 
creased a better sampling is achieved, and the anisotropy 
vanishes. We find that a 216 x 214 A^ simulation cell, 
containing 17200 atoms, is required to obtain a converged 
value K = 1015 ± 120 Wm~^K~^. Our estimate of k is 
lower than the values reported in recent works, in which 
smaller systems were simulated^"—. The inset in Fig. 
[1] displays the normalized heat flux autocorrelation func- 
tions (HFACF) for several simulations with different cell 
sizes, showing that in all cases the time decay is faster 
than 1/t, which guarantees convergence of Eq. ([T]). We 
can conclude that k of unstrained graphene at finite tem- 
perature is finite and converges with size, confirming the 
prediction of former ab initio LD calculationsi^. 

Classical calculations of k far below the Debye temper- 
ature (8d ^ 2000A' for graphene) may yield large dif- 
ferences with respect to calculations taking into account 
the proper quantum statistics for phonons. Two effects 
contribute to such differences, yet in opposite directions: 
classical calculations give shorter phonon lifetimes than 
quantum calculations, but classical phonon heat capac- 
ities are always larger than quantum ones. LD calcula- 
tions showed that in graphene at room temperature these 
two effects compensate to the point that classical k un- 
derestimates quantum k only by about 10%i^. 

The observed size convergence of k from above, pro- 
vides an insight into the contribution of fiexural phonons 
to thermal transport in graphene. Convergence trends 
can be interpreted by referring to the dispersion rela- 
tions of phonons in graphene, computed in the harmonic 
approximation by diagonalizing the dynamical matrix^i 
(Fig. [2|). Accurate MD calculations of the thermal con- 
ductivity of graphene or carbon nanotubes require a con- 
verged sampling of the low-frequency acoustic fiexural 
(ZA) modest. Whereas the contribution to k of in-plane 
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FIG. 2. (Color online) Dispersion relations of unstrained 
graphene (black) and of graphene under uniaxial strain of 4% 
in the zig-zag direction (red). Dispersion curves are computed 
in the strained armchair direction (F-M) and in the stress-free 
zigzag direction (F-K). Wavevectors are in units of the inverse 
lattice constant 27r/a, rescaled according to applied strain. 
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FIG. 3. (Color online) Average thermal conductivity of 
graphene in the presence of flexural vibrational modes (3D), 
and in the absence of such modes (2D), calculated from the 
Green-Kubo relation. The dashed horizontal lines indicate 
the uncertainty for the case of 3D graphene. The inset shows 
the normalized heat flux autocorrelation functions, which de- 
cay like Ijt for the 2D case, yielding log(t) divergence of k. 



acoustic modes converges relatively fast, good sampling 
of the ZA modes is achieved only for large simulation cells 
because of their quadratic dispersion relation near the F 
point. ZA modes are expected to provide a significant 
contribution to heat transport, and have been identified 
as the majority heat carriersi^. However, close to the F 
point their group velocity vanishes and their main role in 
thermal transport is to scatter other heat carriers. Our 
convergence trends indeed show that the overall effect of 
low-frequency ZA modes is to lower the in-plane thermal 
conductivity of graphene. In fact, performing simulations 
on '2D graphene', i.e., a graphene sheet in which atoms 
move only in plane, we observe logarithmic divergence of 
K(t), as shown in Fig. O in accordance with theoretical 
and numerical studies on 2D model systema^^i"— . Our re- 
sults demonstrate the dual role of ZA modes, which is to 
provide an important reservoir of heat carriers, as well as 
the main scattering channel that prevents the divergence 
of 



B. Thermal conductivity of strained graphene 

Strain affects the vibrational properties of materials, as 
it modifies phonon dispersion relations. Speed of sound, 
frequency range, scattering rates, and therefore thermal 
conductivity, are all altered. We apply uniaxial tensile 
strain along the armchair direction and relax the sim- 
ulation cell to achieve zero stress in the perpendicular 
(zigzag) direction. The dispersion relations of strained 
graphene (strain e = 4%) are compared to the unstrained 
ones in Fig. [51 In the low-frequency range the most sig- 
nificant changes is the linearization of the dispersion re- 
lation of the ZA mode along the strain axis. In the di- 
rection perpendicular to the strain axis the ZA branch 
remains unchanged. This implies non- vanishing group 
velocity for the ZA modes propagating along the strain 



axis. In addition, in-plane acoustic modes arc slightly 
softened in both directions. In the high-frequency range 
the degeneracy of the zone center optical phonon is bro- 
ken, in accordance with Raman measurements.— 

Fig. |4] shows the thermal conductivity of graphene as a 
function of time, calculated as the argument of the time 
limit in Eq. ((!]) , along the strained direction for strain up 
to 6%. For small strain (e = 1%), the thermal conductiv- 
ity in the strained direction still converges, yet to a larger 
value than in the unstrained case. As strain is increased 
(e > 2% in the figure), the thermal conductivity along 
the strained direction tends to diverge. The same be- 
havior is observed for a larger simulation cell containing 
more than 10^ atoms. The inset in Fig. |4] displays the 
time decay of the respective HFACF. For e < 2%, the 
HFACF decays faster than Xjt. However, at larger strain 
the HFACF decays as such that k diverges as log(t), 
following the standard behavior of transport coefficients 
in 2D systems^!. Meanwhile, the thermal conductivity 
in the stress- free direction does not diverge^. In fact, 
the thermal conductivity perpendicular to the strain di- 
rection is slightly reduced, due to a mild softening of the 
in-plane acoustic modes as shown in Fig. [2] 

It is important to point out that a logarithmic diver- 
gence of K with size in finite 2D model systems under 
stationary non-equilibrium conditions implies a logarith- 
mic divergence of k as a function of time in periodic sys- 
tems at equilibrium and vice versed"—. In other words, 
a 1/t decay of the heat flux autocorrelation function in a 
periodic system at equilibrium (in the absence of a tem- 
perature gradient), implies a logarithmic divergence of 
K with system size under non-equilibrium conditions (in 
the presence of a finite temperature gradient). Therefore, 
our predictions can (in principle) be probed experimen- 
tally by measuring the size dependence of the thermal 
conductivity in strained graphene samples. 
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FIG. 4. (Color online) Thermal conductivity of strained 
graphene at (a) 300 K and (b) 800 K, calculated from Eq. 
IJ]). The dashed horizontal lines indicate the average k for 
unstrained graphene. The inset shows the time decay of the 
normalized HFACF. For high strain, the HFACF decays as 
1/t such that K diverges as log(t) as observed in the main 
panel. 



Ab initio LD calculations predict divergence of n in 
isotropically strained graphene for any amount of ap- 
plied straiiii^. MD simulations of isotropically strained 
graphene at 300 K suggest that k diverges already for 
e = 1%, confirming the predictions from LD^. In con- 
trast, uniaxial strain and finite temperature limit diver- 
gence to relatively large strain e > 2%, whereas at lower 
strain k remains finite. Simulations at T — 800 K con- 
firm the divergence of n along the direction of strain 
persists for e > 2%, evidencing no significant difference 
with respect to the trends observed at room tempera- 
ture. In this higher temperature regime, even though 
far below the Debye temperature of graphene, classical 
phonon populations approach quantum populations, and 
quantum effects are mitigated. Since phonon lifetimes 
computed with classical statistics are usually underesti- 
mated with respect to those obtained with the correct 
quantum statisticsi^, we can safely argue that the diver- 
gent nature of k in strained graphene is not an artifact 
of classical MD. 

From both analytical models and LD calculations, it 
appears that the cause for the divergence of k lies in the 
linearization of long- wavelength Z A modes (Fig. ^ . Such 
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FIG. 5. (Color online) Lifetimes of the out-of-plane (ZA) 
phonon modes of suspended graphene model (4400 atoms) 
with 0, 2 and 4% uniaxial strain along the armchair direction, 
at 300 K. For strained samples only the modes propagating 
along the strain direction (F-M) are shown. Lifetimes at small 
wavevectors q diverge for the 4% strained model. The 1/q 
decay is shown (dashed line) for reference. 



alteration of the dispersion relations affects phonon pop- 
ulations, which are probed by computing the vibrational 
density of states (VDOS) of strained and unstrained sam- 
ples. Our calculations indeed show a depletion of the 
VDOS for the ZA modes in the frequency range from 0.1 
THz up to « 15 THz**". These changes have dramatic 
effects on phonon lifetimes (Fig. [5|), which diverge for 
large strain and lead to divergence of the thermal con- 
ductivity. Using the solution of the linearized Boltzmann 
transport equation in the relaxation time approximation, 
one can express k for a periodic system of finite size as 
the sum of the contributions of each phonon mode, as 
K = J^f"'' Civfri. Here q is the specific heat of mode i, 
Vi its group velocity, and Ti its lifetime. Since Vi and Ci 
are always finite, divergence of k implies divergence of Ti 
for some of the modes. Phonon lifetimes are computed 
here as the decay time of the autocorrelation function of 
the energy of the normal modes in microcanonical MD 
simulations^. We indeed observe the presence of ZA 
modes with slowly decaying correlation functions when e 
is larger than 2%, which imply diverging r. The trends 
of T(q) (Fig. [5]) permit the extrapolation of our results 
to extended systems, for which k is expressed in integral 
form: 



K(x I dq c(q)v (q)r(q) 



(2) 



where the integral is taken over the two-dimensional Bril- 
louin zone (BZ). The limit for g — > determines whether 
K diverges. In the classical case c(q) = ks/V is a con- 
stant, and in general t oc l/w". The group velocity of 
the ZA modes along x is v^^ = duj/dqx- For unstrained 
graphene v^^ oc Qx , whereas when strain is applied along 
X, tends to a constant value. In unstrained graphene 
K would diverge for a > 2, while when strain is applied 
K diverges for a > 1. Our calculations show that a in- 
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FIG. 6. (Color online) Temperature dependence of the 
thermal conductivity of unstrained graphene with different 
amounts of ^'^C. Solid lines are fitted to the data in order to 
extract the temperature dependence, as reported in the plot. 



creases with strain, reaching a ^ 1 for e = 2%, which is 
consistent with the observed threshold for k divergence. 



C. Thermal conductivity of isotopically modified 
graphene 

So far we have presented results for isotopically 
pure graphene (100% ^^C). Given the demonstration of 
graphene growth with customized isotopic composition^, 
it is also worth investigating the combined effect of strain 
and controlled isotopic composition on the thermal con- 
ductivity of graphene, checking whether its divergence 
can be suppressed by mass disorder. We consider pure 
^^C, natural composition (1.1% ^^C), and ^'^C enriched 
graphene models (10%, 50% and 99.2%). In absence 
of defects the lattice thermal conductivity is limited by 
phonon-phonon scattering^, therefore k oc l/T. Mean- 
while, the scattering of phonons by defects is tempera- 
ture independcnlj^, so that k{T) trends are modified. In 
Fig. |6] we indeed observe k ^ l/T for isotopically pure 
graphene, natural graphene and 99.2% ^^C isotopically 
enriched graphene, indicating that for the natural iso- 
topic composition the effect of mass disorder is almost 
negligible. However, as the amount of ^'^C increases to 
10% and 50%, k decreases more slowly with T indicat- 
ing that mass disorder becomes the primary source of 
phonon scattering. This aspect could be exploited in 
thermal management devices that operate over wide tem- 
perature ranges, where it would be undesirable to have 
large variations in k with T. The ratio between n for iso- 
topically pure and 50% ^^C-enriched graphene is about 2, 
and agrees well with recent experimental measuremcntsS.. 

Even though large variations of k as a function of iso- 
topic composition are found, when uniaxial tensile strain 
is applied the general behavior of k for isotopically mod- 
ified graphene is not qualitatively different from the iso- 
topically pure case, k along the strained direction in- 
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FIG. 7. (Color online) Thermal conductivity of unstrained 
and strained isotopically modified graphene at 300 K, calcu- 
lated from Eq. ((!]), with a 50% concentration of ^^C. The 
dashed horizontal lines indicate the average k for the un- 
strained case. The inset shows the normalized heat flux auto- 
correlation functions, which decay like 1/t for e > 2%, yield- 
ing log(t) divergence of k. 



creases for low strain and diverges for e > 2% even in 
samples with the highest isotopic disorder (50% ^^C), as 
shown in Fig. [T] The inset shows that above the strain 
threshold the HFACF decay as 1/t and thus k ^ log(t). 
The ratio between phonon populations shows that in iso- 
topically enriched graphene, as in the isotopically pure 
case, tensile strain induces similar reductions to the pop- 
ulation of ZA modest. Therefore, mass disorder is not 
sufficient to suppress the divergence of the lifetimes of 
low-frequency ZA modes. In fact, these modes have a 
wavelength of several tens of nm, thus even in the pres- 
ence of isotopic mass disorder they propagate as in a 
continuous medium, and are not significantly affected by 
scattering centers at the atomic scale. It is worth not- 
ing that mass disorder can not suppress divergence in 
non-linear models as well^ii^. 



IV. CONCLUSIONS 

In conclusion, we have shown that heat transport in 
suspended graphene is controlled by ZA modes, which 
contribute the essential scattering channels to limit the 
thermal conductivity in unstrained samples. In fact, in 
absence of ZA modes k would diverge with size, as in 
ideal 2D models. Uniaxial tensile strain reduces the pop- 
ulation of ZA phonon modes at low frequency, makes 
their zone center group velocity finite and increases their 
lifetime, thus causing divergence of the thermal conduc- 
tivity in the strained direction. We then predict that k 
of strained samples would diverge logarithmically with 
the size of the samples. It is important to point out 
that our predictions based on computer simulations are 
accessible to experiments. In an experimental setup, k 
will always be finite and limited by boundary and defect 
scattering. Nonetheless, by performing measurements in 



6 



strained samples of increasing size it should be possible 
to see a logarithmic dependence of function of 

size. The amount of strain required to observe diver- 
gence, £ > 2%, might also be within reach of current 
experimental techniques. We also predict that the pres- 
ence of isotopic mass disorder does not suppress the di- 
vergence of K that should be then expected to occur also 
in samples with natural composition, which are easier to 
grow than isotopically pure ones. 
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